## Code for "Dynamic Parties and Social Turnout"
## James Fowler
## 11/7/03

rm(list=ls(all=TRUE))

DPAST<-function(n=32,p=c(-1,1),cost=0.1,vMem=0.0,prefCor=0.0,
  partiesFixed=FALSE,turnoutFixed=FALSE,B=1,nPeriods=100,nSims=1,
  graphs=FALSE,fileName=NULL) {

## key parameters
## n<-32  ## length, width of grid
## p<-c(-1,1)  ## party preferences
## cost<-0.1 ## cost of turnout
## vMem<-0.0 ## voter memory weights past vs. present observations
## prefCor<-0.75  ## degree of correlation of voter preferences
## partiesFixed=FALSE  ## choose whether or not to fix parties
## turnoutFixed=FALSE  ## choose whether or not to fix turnout
## B<-1  ## voter polarization (var of voter distribution)
## nPeriods<-100  ## number of periods
## nSims<-1  ## number of repetitions of simulation

for (sims in 1:nSims) {

## intialize time series and other variables
N<-n^2  ## number of voters
x<-array(0,dim=c(nPeriods,2)) ## platforms
## prior mean, variance, turnout, vote share, median voter
mu<-bt<-turnout<-s<-mv<-rep(0,nPeriods) 

## initialize count arrays
nt<-npt<-npa<-array(0,dim=c(n,n))

## set intial turnout unless fixed
iTurnout<-0.5 ## initial turnout
if (turnoutFixed) iTurnout=1.0

## set initial party platforms
x[1,]<-p

## initialize voter preferences as combination of sorted and 
## unsorted preferences depending on preference correlation
sorted<-rbinom(N,1,prefCor)
pref<-array(sorted*sort(rnorm(N,0,sqrt(B)))+
            (1-sorted)*rnorm(N,0,sqrt(B)),dim=c(n,n))

## FOCs for left and right party
deuL <- function(x,pL,mu,bt) -(exp(-((sum(x)-2*mu)^2)/(8*bt))*
  (x[1]-x[2])*(-2*pL+sum(x)))/(2*sqrt(2*bt*pi))+
  2*(pL-x[1])*pnorm(mean(x), mean=mu, sd=bt^(1/2))
deuR <- function(x,pR,mu,bt) -(exp(-((sum(x)-2*mu)^2)/(8*bt))*
  (x[1]-x[2])*(-2*pR+sum(x)))/(2*sqrt(2*bt*pi))+
  2*(pR-x[2])*(1-pnorm(mean(x), mean=mu, sd=bt^(1/2)))

## function for collecting neighborhood stats
neighbors<- function(mat) {
  mat<-rbind(rep(0,(n+2)),cbind(rep(0,n),mat,rep(0,n)),rep(0,(n+2)))
  mat[1:n,1:n]+mat[1:n,2:(n+1)]+mat[1:n,3:(n+2)]+
    mat[2:(n+1),1:n]+mat[2:(n+1),2:(n+1)]+mat[2:(n+1),3:(n+2)]+
    mat[3:(n+2),1:n]+mat[3:(n+2),2:(n+1)]+mat[3:(n+2),3:(n+2)]
}

## index array of number of neighbors in each neighborhood
nn<-neighbors(array(1,dim=c(n,n)))

## INITIAL ELECTION

  ## intialize turnout actions
  action<-array(rbinom(N,1,iTurnout)==1,dim=c(n,n))

  ## voters choose parties (flip coin if equidistant)
  choice<- (-(x[1,2]-pref)^2 > -(x[1,1]-pref)^2) +
           (-(x[1,2]-pref)^2 == -(x[1,1]-pref)^2)*rbinom(N,1,0.5)

  ## determine turnout, vote share, winner (flip coin if tie)
  turnout[1]<-mean(action)
  if (turnout[1]>0) {
    s[1]<-mean(action*choice)/turnout[1]
  } else {
    s[1]<-0.5  ## ignorance result when no one votes
  }
  s[1]<-max(min(s[1],(N-1)/N),1/N)  ## restrict s[t] to avoid inf in m

  winner<-(s[1]>0.5)+(s[1]==0.5)*rbinom(1,1,0.5)

  ## voter payoffs
  payoff<- -(x[1,winner+1]-pref)^2-action*cost

## ELECTION CYCLE
for (t in 2:nPeriods) {

  ## PARTIES

  if (partiesFixed) {
     x[t,]<-x[t-1,]

  } else {
    ## infer median voter from previous vote share and platforms
    m <- qnorm(s[t-1],mean(x[t-1,]),sqrt(B))

    ## Bayesian update beliefs
    mu[t]<-((t-2)*mu[t-1]+m)/(t-1)  ## mean
    bt[t]<-((t-2)*bt[t-1]+(m-mu[t-1])^2)/(t-1)  ## variance

    ## optimize by minimizing loss function
    peq<-optim(c(mu[t],mu[t]+.01),function(z)
      abs(deuL(z,p[1],mu[t],max(bt[t],.000001))) + 
      abs(deuR(z,p[2],mu[t],max(bt[t],.000001))) +
      if(z[1]<p[1]||z[2]>p[2]) .01 else 0)

    ## set eq platforms
    x[t,]<-peq$par
  }

  ## VOTERS

  if (!turnoutFixed) {

    ## count number of voters in each neighborhood
    nt<-vMem*nt+(1-vMem)*neighbors(action)

    ## sum payoffs of voters and abstainers in each neighborhood
    npt<-vMem*npt+(1-vMem)*neighbors(payoff*action)
    npa<-vMem*npa+(1-vMem)*neighbors(payoff*(1-action))

    ## turn out if ave payoff of voting neighbors is higher than 
    ## ave payoff of abstaining neighbors
    action<-((npt/nt)>(npa/(nn-nt))) 
    action<-action&(nt>0)  ## turn out only if someone else turns out
    action<-action|(nn-nt==0)  ## turn out if no one else abstains
  }

  ## voters choose parties (flip coin if equidistant)
  choice<- (-(x[t,2]-pref)^2 > -(x[t,1]-pref)^2) +
           (-(x[t,2]-pref)^2 == -(x[t,1]-pref)^2)*rbinom(N,1,0.5)

  ## determine percent turnout, vote share, winner (flip coin if tie)
  turnout[t]<-mean(action)
  if (turnout[t]>0) {
    s[t]<-mean(action*choice)/turnout[t]
  } else {
    s[t]<-0.5  ## ignorance result when no one votes
  }
  s[t]<-max(min(s[t],(N-1)/N),1/N)  ## restrict s[t] to avoid inf in m
  winner<-(s[t]>0.5)+(s[t]==0.5)*rbinom(1,1,0.5)

  ## median voter 
  mv[t]<-median(pref[which(action)])

  ## voter payoffs
  payoff<- -(x[t,winner+1]-pref)^2-action*cost
}

## median preference
medPref<-median(pref)

if (graphs) {
  par(mfrow=c(2,3)) ## clear graphs
  plot(x[,1], 
type="l",ylim=c(min(c(x[2:nPeriods,],mv)),max(c(x[2:nPeriods,],mv))))
  lines(x[,2], type="l", lty=2)
  plot(mv, 
type="l",ylim=c(min(c(x[2:nPeriods,],mv)),max(c(x[2:nPeriods,],mv))))
  lines(rep(median(pref),nPeriods), type="l", lty=2)
  plot(mu, 
type="l",ylim=c(min(c(x[2:nPeriods,],mv)),max(c(x[2:nPeriods,],mv))))
  plot(turnout, type="l")
  plot(s, type="l")
  lines(rep(0.5,nPeriods), type="l", lty=2)
  plot(bt, type="l")
}

## turnout correlation
turnoutCor<-cor(c(action[2:n,2:n],action[2:n,1:n],action[2:n,1:(n-1)],
  action[1:n,2:n]),c(action[1:(n-1),1:(n-1)],action[1:(n-1),1:n],
  action[1:(n-1),2:n],action[1:n,1:(n-1)]))

## preference correlation
prefCorr<-cor(c(pref[2:n,2:n],pref[2:n,1:n],pref[2:n,1:n-1],pref[1:n,2:n]),
  c(pref[1:(n-1),1:(n-1)],pref[1:(n-1),1:n],pref[1:(n-1),2:n],
  pref[1:n,1:(n-1)]))

## interior exterior turnout
intTurnout<-mean(action[2:(n-1),2:(n-1)])
extTurnout<-mean(c(action[1,],action[,1],action[n,],action[,n]))

## write data if file name provided
if (!is.null(fileName)) 
write.table(t(c(n,p,cost,vMem,prefCorr,partiesFixed,
    
turnoutFixed,B,nPeriods,turnoutCor,prefCor,intTurnout,extTurnout,medPref,
    x,turnout,s,mv,mu,bt)),file=fileName,append=TRUE,row.names=FALSE, 
col.names=FALSE)


} ## end sim
} ## end function